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Chapter  I 


INTRODUCTION 


Many  problems  in  continuum  mechanics  and  mechanical 
system  design  involve  elastic  bodies  that  come  into  con- 
tact with  each  other,  under  applied  loads.  Such  prob- 
lems, called  contact  problems  of  elasticity,  are  non- 
classical  in  the  sense  that  one  does  not  initially  know 
the  contacting  region.  Considerable  research  has  been 
pursued  in  recent  years  to  develop  constructive  methods 
of  determining  the  contact  region  and  contact  stress 

K 

distribution  [l-4]. 

The  geometry  of  two  bodies  that  will  come  into  con- 
tact under  applied  load  is  shown  schematically  in  Fig- 
ure 1(a).  A numerical  method  of  solution,  by  quadratic 
programming  [5]  has  been  developed,  in  which  one  selects 
a set  of  potential  contact  points,  as  indicated  in  Fig- 
ure 1(a).  The  selection  of  potential  contact  points  is 
explained  in  detail  in  [6],  The  initial  gap  between  the 
ith  potential  contact  points  is  denoted  by  d^.  The  rig- 
id body  displacement  coordinates  of  Body  One  are  indi- 
cated in  Figure  1(a)  and  are  represented  by  a general- 
ized coordinate  vector  q,  of  dimension  one  to  six. 

Once  the  potential  contact  points  are  defined  on  the  two 
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bodies,  contact  stress  may  be  replaced  by  equivalent  point 
loads,  or  contact  forces,  denoted  as  in  Figure  1(b). 

When  the  bodies  come  into  contact,  a force  distribu- 
tion over  the  surfaces  arises  and  high  contact  stresses  may 
occur  over  some  parts  of  the  contact  region.  This  is  unde- 
sirable, since  plastic  deformation  of  the  bodies  may  occur, 
or  high  normal  forces  may  lead  to  wear  of  machine  parts  that 
move  relative  to  each  other.  This  can  be  quite  important  in 
precision  machines.  The  objective  of  this  work  is  to  devel- 
op a technique  for  adjusting  the  contour  of  one  or  both  of 
the  bodies  in  order  to  achieve  a minimum  peak  contact  stress 
between  the  bodies,  under  a given  load. 

A related  problem  was  solved  by  Conry  [6],  in  which  the 
design  objective  was  to  select  the  contour  to  achieve  a con- 
stant stress  over  the  contact  arc.  A linear  programming 
method  was  presented  in  [7,8]  to  implement  this  design  ob- 
jective. In  general,  one  cannot  achieve  a constant  stress 
distribution  over  the  contact  arc  when  kinematic  constraints 
define  limits  on  modification  of  the  contours  of  the  bodies. 
This  will  be  particularly  critical  in  precision  cams  and 
fasteners  that  require  precise  location  of  parts.  The  de- 
sign objective  selected  in  this  work  is  to  minimize  the  peak 
contact  stress,  subject  to  constraints  on  the  extent  of 
modification  of  surface  contours. 


Body  2 


1(a). 


1(b). 


Figure  1 . 


Geometry  of  Contacting  Bodies. 
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Chapter  II 


THE  UNBONDED  CONTACT  PROBLEM  FOR 
ELASTIC  BODIES  IN  CONTACT 


2 . 1 Analytical  Formulation  of  the  Problem 

The  elastic  contact  problem  formulation  of  [,6^  will 
be  employed  as  the  model  of  the  elastic  contact  problem. 
Only  a brief  summary  of  this  analytical  formulation  will 
be  presented  here.  The  gap  variable  £,  due  to  a stress 
distribution  Sf  is  given  by  the  equation 

£ = BS  + Aq  + a + b (2.1.1) 

where  B is  a matrix  formed  from  influence  coefficients 
of  the  two  bodies,  A is  an  affine  transformation  account- 
ing for  rigid  body  displacement  of  Body  1,  a is  a con- 
stant vector  that  depends  on  the  externally  applied  loads 
to  the  two  bodies  and  the  initial  gap  between  them  at 
the  potential  contact  region,  and  b is  the  vector  of 
contour  modification  to  be  selected.  Formulae  for  B, 

A,  and  a are  given  in  Appendix  A. 

Equilibrium  of  Body  l is  obtained  through  direct 
application  of  the  principle  of  virtual  work,  which  re- 
sults in  the  linear  equation 
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ATS  = C (2.1.2) 

where  the  column  vector  C depends  on  the  externally  ap- 
plied load,  and  is  given  in  Appendix  A. 

Compatibility  conditions  between  the  two  contacting 
bodies  require  that  the  product  of  each  gap  variable  and 
contact  stress  be  zero.  Analytically,  this  is 

€iSi  = 0 , i = 1,2, , n (2.1.3) 

where  n is  the  number  of  potential  contact  points. 

This  condition  may  be  interpreted  as  stating  that  either 
the  gap  or  the  contact  stress  must  be  zero  at  each  point 
on  the  contacting  bodies.  Further,  the  gap  and  contact 
stress  must  be  non-negative.  Analytically  these  condi- 
tions are 

Sx  ^ 0,  % 0,  i = 1,2,  , n (2.1.4) 

As  is  shown  in  [6],  Equations  (2.1.l)-(2.l.4)  are 
the  Kuhn-Tucker  necessary  conditions  for  solution  of  a 
convex  quadratic  programming  problem,  which  may  be 
stated  as 

minimize  ATS  + i gT^ 

subject  to  constraints  _ c 

S £ 0 
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This  problem  may  be  solved  through  direct  application  of 
the  simplex  technique  of  quadratic  programming  [5],  It 
is  also  shown  in  [5]  that  if  a solution  exists,  it  is 
unique  and  may  be  reliably  determined,  through  applica- 
tion of  a simplex  technique  for  quadratic  programming 
[6]. 

The  cost  function  to  be  minimized  in  this  design 
problem  is  the  maximum  contact  stress, 

J = Max  (2.1.6) 

Since  this  maximum  value  function  is  difficult  to 
treat  analytically,  it  may  be  replaced  by  defining  an 
upper  bound  bfi  + on  the  contact  stress,  through  the 
set  of  inequalities 

$i  = Si  - bn  + 1*  0 <2-1'7) 

The  cost  function  J of  equation  (2.1.6)  may  now  be  re- 
placed by  an  equivalent  problem  of  minimizing  the  upper 
bound  on  contact  stress 

J = b . (2.1.8) 

n + l 


subject  to  constraints  (2.1.7). 


Finally,  it  is  often  required  that  both  upper  and 
lower  bounds  be  placed  on  modifications  of  the  surface 
contours  of  the  two  bodies.  Analytically,  these  condi- 
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tions  may  be  stated  as 

^>n  + ^ = b°  0 , i = 1, , n (2.1.9) 

and 

^2n  + i = b^  - b*  ^ 0 , i = 1, , n (2.1.10) 

where  b°  and  b1  are  the  lower  and  upper  bounds  on  contour 
modification,  respectively. 

The  design  problem  may  now  be  stated  as  follows: 
Determine  the  design  variable  vector  b and  the  upper 
bound  bfi  + 1 to  minimize  J,  subject  to  constraints 
(2.1.1)  - (2.1.4)  , (2.1.7),  (2.1.9)  and  (2.1.10). 

2 . 2 Sequential  Linear  Programming 
Solution  Technique 

With  an  initial  estimate  of  the  contour  (normally 
b = 0),  one  may  solve  the  contact-analysis  problem  of 
equations  (2.1.1)  - (2.1.4)  to  obtain  the  contact  arc 
and  contact  stresses.  Using  this  solution,  define  the 


index  sets 


A 

I 


i 


8 


of  points  in  contact  and 


i 


j 


of  points  not  in  contact,  where  NNC  is  the  number  of 
points  not  in  contact.  The  stresses  at  points  in  con- 

A 

tact  will  be  denoted  as  the  NC  dimensional  vector  S. 

A 

Define  B as  the  matrix  B with  rows  and  columns 

A 

corresponding  to  points  not  in  contact  removed,  A as 
the  matrix  A with  rows  corresponding  to  the  points  not 

A A 

m contact  removed,  and  the  vectors  a and  b as  the 
vectors  a and  b with  components  corresponding  to 
points  not  in  contact  removed. 

One  may  now  state  the  set  of  conditions  that  must 
be  satisfied  if  the  contact  arc  is  not  changed;  namely 

A aA  A a a 

S = BS  + Aq  + a + b = 0 (2.2.1) 


The  equilibrium  equations  (2.1.2)  must  still  be  satis- 
fied with  contact  loads  applied  only  on  the  fixed  con- 
tact region,  so  one  obtains 
a m a 

A s = C (2.2.2) 


The  objective  now  is  to  determine  the  reduced  design 

A 

variable  b so  that  peak  contact  stresses  are  reduced 
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as  much  as  possible,  keeping  the  same  contact  region  and 
satisfying  conditions  (2.1.7)  - 2.1.10). 

The  gap  variable  outside  the  prescribed  contact  re- 
gion must  remain  non-negative,  so 


where  B is  the  matrix  B with  columns  corresponding 
to  points  not  in  contact  and  rows  corresponding  to  the 


points  in  contact  removed,  A is  the  matrix  A with 
rows  corresponding  to  points  in  contact  removed,  and 
vectors  a and  b are  the  vectors  a and  b with 
rows  corresponding  to  points  in  contact  removed.  Here, 

AS 

the  reduced  vector  b represents  the  estimates  of  the 
contour  design  variable  from  the  preceding  iteration. 
Finally,  it  is  necessary  to  impose  the  condition. 


The  design  modification  problem  is  now  to  deter- 


Af  fSA.  A/  r*  A*  ^ 

6 = BS+Aq+a  + b>0 


(2.2.3) 


(2.2.4) 


mine  the  reduced  vector  b and  b 


n + 1 


to  minimize 


J = b 


n + 1 


subject  tot 


(2.2.5) 
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§ ^ b 


b°  - b 


A 

b - b 


n + 

i < 


i 

0 


aa  a 

BS  + Aq  + 

A T A 
A 1 S = C 


A 

a 


+ b = 0 


(.  <*/A  A/  V /V 

8 = BS  + Aq  + a + b 


A 

s 


£ 


0 

0 


(2.2.5,  con ' d . ) 


This  formulation  for  determining  the  modified  contour  var- 
iable is  simply  a linear  programming  problem  in  the  vari- 

A A 

ables  s,  q,  and  b which  can  be  solved  with  standard  lin- 
ear programming  codes. 

Following  each  solution  of  the  linear  programming 
problem,  the  contact  stresses  and  gap  variables  must  be 
evaluated  to  determine  the  contact  surface  to  be  employed 
in  the  next  iteration.  Any  potential  contact  point  for 
which  Sj  = 0 is  included  in  the  contact  surface  for  the 
next  iteration.  At  any  points  in  the  preceding  potential 

A 

contact  surface  for  which  Sj  = 0,  it  is  presumed  that 
separation  will  occur  in  the  subsequent  iteration,  so 

A A 

these  points  are  deleted  from  the  vectors  8 and  S . 
These  modification  rules  form  the  basis  for  definition 
of  the  next  linear  programming  problem  to  be  solved  for 
the  optimum  contour  on  an  adjusted  contact  arc. 

This  process  is  continued  until  no  new  points 
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come  into  contact  and  no  points  that  were  previously  in 
contact  are  released. 

Each  linear  programming  solution  will  reduce  the 
value  of  the  peak  stress,  until  the  process  stops. 
Analytically,  this  condition  can  be  written  as  follows 


b 

b 


(o) 
n + 1 

£ 

b(1>  , 

n + 1 

(1) 

n + 1 

£ 

»(2i , 

n + 1 

(k-l) 


>(*> 

n + 1 


(2.2.6) 


where  the  superscript  on  br 
number.  Combining  (2.2.6) 


>(°) 
n + 1 


b{l)  > 
n + 1 x 


denotes  the  iteration 


one  obtains 


* n + 1 (2.2.7) 


This  gives  a sequence  of  non-increasing  real  numbers, 
called  a monotone  sequence,  that  is  bounded  below  by 
zero.  Any  non-increasing  sequence  that  is  bounded  be- 
low is  convergent  [9j,  so  the  sequence  of  solutions  of 
linear  programming  problems  (2.2.5)  must  converge. 

It  was  observed  from  preliminary  calculations  that 
as  the  value  of  U is  increased,  peak  stress  decreases, 
but  a stage  is  reached  when  stress  at  all  the  points  in 
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the  contact  region  is  the  same.  The  peak  stress  obtained 
in  this  case  is  a local  minimum  and  no  contour  design  vari- 
able reaches  its  allowable  limit.  This  situation  can  be 
avoided  by  adding  a penalty  function  to  the  cost  function 
that  is  intended  to  broaden  the  contact  region.  The 
penalty  function  used  here  is  as  follows 


J = b 


n + 1 


„ NNC 
*1  * 


J = 1 


e. 

j 


(2.2.8) 


where  is  a small  constant  greater  than  zero.  The  linear 
programming  with  the  augmented  cost  function  (2.2.8)  finds 
an  absolute  minimum  of  the  peak  contact  stress. 

The  process  described  above  is  summarized  in  the 
following  algorithms 


Step  1.  Estimate  the  design  variables  (normally  b^°^ 
= 0)  . 

Step  2.  Solve  for  S and  o by  quadratic  programming, 
using  the  necessary  conditions  (2.1.1)- 
(2.1.4) . 

A a/ 

Step  3.  Form  index  sets  I and  I for  points  of  con- 
tact (S^  ^ 0,  £i  = 0)  and  for  points  not  in 
contact  ( ^ 0,  S^  = 0),  respectively. 

A ^ A 

Step  4.  Construct  B,  A,  a,  and  b for  lei  and  A, 


B,  a,  and  b for  jel. 
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Step  5. 

Step  6. 
Step  7. 


Step  8. 


Solve  the  linear  programming  problem  (2.2.5) 

A A . 

for  S,  q,  and  b.  (The  cost  function  may  be 
changed  as  given  in  (2.2.8)  if  necessary). 

A/ 

Evaluate  £ from  (2.2.3). 

From  new  I to  include  points  for  which  Cj 
= 0 and  deleting  points  for  which  Sj  = 0 
and  form  I of  points  for  which  Oj  > 0 and 

Sj  s 0. 

A 

If  I is  unchanged,  terminate;  otherwise  re- 
turn to  Step  4. 


Chapter  III 


CONTACT  SURFACE  DESIGN  PROBLEMS  FOR 
BEAMS  ON  ELASTIC  FOUNDATIONS 


Examples  are  presented  to  demonstrate  the  tech- 
nique developed  in  Chapter  II.  Two  cases  of  the  design 
problem  with  a beam  on  an  elastic  foundation  are  con- 
sidered! (l)  with  no  initial  gap;  and  (2)  with  an  ini- 
tial gap. 

A point  load  is  applied  to  the  center  of  the  beam, 
as  shown  in  Figure  2(b).  Only  vertical  displacement  of 
the  center  of  the  beam  is  chosen  to  represent  its  rigid 
body  displacement.  Formulation  of  the  matrices  B,  A,  a 
b,  and  C for  this  problem  is  explained  in  detail  in  [2] 
A quadratic  programming  problem  is  solved  to  obtain  the 
contact  arc  as  input  to  the  linear  programming  problem. 

The  design  variables  b.  are  limited  by  the  con- 

Xj2  1 

straint  b U (l  - — ~ ),  where  x.  is  the  horizontal 

1 Lz  1 

distance  of  the  ith  potential  contact  point  from  the 
center  of  the  beam,  L is  half  the  potential  contact 
length  (as  only  half  of  the  beam  is  considered  due  to 
symmetry),  and  U is  a constant  greater  than  zero. 
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Figure  2. 


Potential  Contact  Points  for  the 
Beam  on  Elastic  Foundation. 
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3 . 1 Beam  on  Elastic  Foundation  with  No  Initial  Gap 


The  following  material  properties  are  used:  Young's 

modulus  of  the  beam  material  = 340,000  psi ; diagonal  ele- 
ment of  the  flexibility  matrix  for  the  elastic  foundation 
=0.00005  in. /lb. ; height  of  the  beam  = 0.5  in.;  and  width 
of  the  beam  = 1.0  in. 

The  quadratic  programming  problem  was  solved  for  25 
potential  contact  points  on  half  of  the  beam,  with  an  in- 
terval size  of  0.06  in.  Potential  contact  length,  as  shown 
in  Figure  3(a)  is  taken  as  2.88  in.  The  contact  length  ob- 
tained numerically  was  1.38  in.,  with  loads  of  1000  lbs. 
and  2000  lbs.  A quadratic  programming  code  called  ZORILLA 
L 1 0 ’ was  used  to  solve  this  quadratic  programming  problem. 
Results  are  presented  in  Figure  3(b). 

At  this  stage,  one  has  the  contact  arc  to  be  employed 
in  the  linear  programming  problem  of  (2.2.5).  Equations 
(2.2.5)  can  now  be  written  as: 


J = b 

A 

S - b 


n + 1 

n + 1 ^ 0 


b 

A 


A A 

BS 

A r 

A ' 


A 

s 


Aq  + b = 
= C 


A 

-a 


/V  A a/  ^ a/ 

BS  + Aq  > -a-b 

A 

s 


(3.1.1) 


? 0 


Contact  Stress  (lbs./sq. 


¥ 


^ Contact  Length  L (in.) 

(b) 

Figure  3.  Contact  Stress  vs.  Contact  Length  from  Quadratic 

Programming  Solution. 
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The  linear  programming  subroutine  used  to  solve  (3.1.1) 

is  explained  in  [ll]  and  is  summarized  in  Appendix  B. 

This  subroutine  requires  that  the  right  hand  side  values 

of  constraints  must  be  given  as  a separate  vector,  all 

2 

elements  of  which  must  be  positive.  If  b°  = -u[l-^~2  ] 

2 u 
and  b1  = u[l  ~ l?  ]>  then  Equations  (3.1.1)  can  be  re- 
written as 


J 

A 

s 

A 

-b 

A 

b 


n + 1 


3n  + 1 < 


AA  A A 

BS  + Aq  + b = 0,  as 

A <p  A 

A 1 S = C 


A 

a 


-BS  - Aq  ^ a r b 


(3.1.2) 


The  linear  programming  subroutine  used  solves  only  max- 
imization problems,  so  the  cost  function  is  replaced  by 

J=  -b  , . Since  q can  be  positive  or  negative, 
n + 1 r 

the  vector  q is  decomposed  into  positive  and  negative 
parts  as  follows: 


q 

+ 

q 

q~ 


q 

q 

o 

q 

o 


- q 


when  q > 0 

q ^ 0 (3.1.3) 

when  q £ 0 

q £ 0 
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Nov  an  LP  tableau  is  formed  as  follows  : 


1 

NC 

NC 

NC 

NC 


m 


NNC 


m m 1 NC  NC 


r — 

— 

— - 

0 ; 0 ; -1  • 0 ; 0 

+ 

1 * * 

, 1 • ' 

q 

. 1 « 1 

H 

O 

M 

1 

O 

O 

i 

0 

• i 

q” 

O 

o 

o 

H 

O 

0 

& 

1 

' . ' 

b 

1 i | | 

n+1 

o 

* 

o 

0 

1 

M 

o 

=■ 

b1 

J- \ j— - 

A 

b 

— 

1 1 1 1 
A | A | i i A 

a i -a  • o ; i ; b 

i ! t i 

0 

, • 

A 

1 ! ! ! a 

S 

0 I 0 ; 0 ; 0 1 A 

' ! 4 I. 

c 

-a  : a j o : o : b 

a+b 

_ _ 

(3.1.4) 


r*  i ' 1 

where  the  vector  NBP  [_  0 j 0 | -1  j 0 j 0 ] of  dimen- 
sion (1,  2m  + 1 + 2 * NC  + NGE)  represents 

the  cost  function  and  NGE  is  the  number  of  constraints 
of  the  form  ^ 0,  and  m is  the  number  of  degrees 


of  freedom.  The  Matrix  B 


defined  as 
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contains  constraint  coefficients.  Its  dimension  is 
(4  * NC  + m + NNC  , 2*m+l+2*NC+  NGE)  and  the 
vector  RQ1  = [ 0 J -b°  • b j o i c i a+b  ],  of  dimen- 
sion (l,  4 * NC  + m + NNC),  is  the  vector  of  right  hand 
side  values  of  the  constraints. 


The  design  problem  was  solved  with  loads  t = 1000 
lbs.  and  t = 2000  lbs.,  with  a range  of  values  of  contour 
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modification  limit  U.  Results  are  presented  in  Figures 
(4  - 9) . All  the  calculations  were  done  on  an  IBM  360/65 
computer.  Table  1 shows  the  comparison  of  results  for 
different  loads  and  values  of  U ranging  from  0.01  to 
0.05.  It  is  observed  that  the  contact  length  for  the  op- 
timum design,  with  a fixed  value  of  U,  varies  with  ap- 
plied load.  However,  when  the  value  of  U is  changed  in 
the  same  proportion  as  the  loads,  it  is  noted  that  the 
contact  length  is  the  same  for  both  cases  (and  peak  con- 
tact stress  are  proportional).  As  the  value  of  U is 
increased,  the  contact  length  appears  to  increase  and  the 
value  of  peak  stress  decreases.  A stage  is  reached  when 
stress  at  each  point  in  the  contact  region  is  the  same, 
which  stops  the  iterative  process.  In  the  case  of  U = 
0.025,  t = 1000  lbs.,  and  U = 0.05,  t = 2000  lbs.,  no 
contour  design  variable's  constraint  reaches  its  tolerance 
limit.  The  cost  function  was  augmented  with  penalty  func- 
tion and  the  linear  programming  problem  was  solved  for  the 
above  case.  The  results  are  presented  in  Figures  (8  - 9). 
It  is  also  observed  that  the  contour  design  variable  con- 
straint at  the  center  point  of  the  contact  region  is  al- 
ways tight. 

Table  1 presents  peak  stress,  computing  time  and  con- 
tact lengths  for  several  numerical  examples.  Solutions  of 
the  linear  programming  problem  converge  in  at  most  8 


Contact  Stress 


Contact  Length  (in.) 


Figure  4.  Contact  Stress  vs.  Contact  Length  with  U = 0.01. 


Tolerance  Limit 


Final  Gap 


> Potential  Contact  Point 


Figure  5.  Gap  Size  vs.  Potential  Contact  Point  with  U = 0.01. 


Contact  Stress 


Figure  6.  Contact  Stress  vs.  Contact  Length 


Potential  Contact  Point 


Figure  7.  Gap  Size  vs.  Potential  Contact  Point 


Contact  Stress  (lbs./sq.  in. 


1200 


800 


400 

t . 


Figure  8.  Contact  Stress  vs.  Contact  Length  with  Augmented  Cost  Function. 


Gap  Size  (in. ) 


Tolerance  Limit 


. Final  Gap 


y Potential  Contact  Point 

Figure  9.  Gap  Size  vs.  Potential  Contact  Point  with  Augmented  Cost  Function. 
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TABLE  1 

COMPARISON  OF  RESULTS  FOR  LOADS 
t = 1000  lbs.  AND  t = 2000  lbs. 
FOR  BEAM  ON  ELASTIC  FOUNDATION 


LOAD 

t 

SR. 

NO. 

U 

PEAK 

STRESS 

lbs./sq. 

COMPUTING 
TIME  (SEC. 
in. 

NO. 

) OF 
ITRNS 

TIME/  CON- 
ITRN  TACT 
(SEC.)  LENGTH 
(IN.) 

1 

.01 

562.9 

38.80 

5 

7.36 

1.86 

2 

.015 

499.2 

48.12 

6 

8.02 

1.98 

1000 

3 

.02 

457.4 

72.08 

8 

9.01 

2.22 

4 

.025 

426.43 

91.27 

9 

10.14 

2.34 

1 

o 

ro 

1125.6 

36.87 

5 

7.374 

1.86 

2 

.03 

998.25 

48.94 

6 

8.157 

1.98 

2000 

3 

.04 

914.8 

73.00 

8 

9.125 

2.22 

4 

.05 

852.87 

93.17 

9 

10.35 

2.34 
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iterations,  with  computing  time  per  iteration  varying  only 
between  7 and  10  seconds.  In  all  cases,  peak  stress  is 
reduced  considerably  and  the  contact  arc  increases.  For 
example,  for  t = 1000  lbs.,  the  peak  stress  of  the  non- 
modified  design  was  1120  lbs. /in.  and  contact  arc  length 
was  1.38  in.  By  adjusting  the  contours  of  two  bodies, 
with  a tolerance  limit  of  U = 0.01  and  t = 1000  lbs.. 

Table  2 shows  the  peak  stress  reduces  to  552.9  lbs./sq.  in. 
and  the  contact  arc  length  becomes  1.86  in. 

3 . 2 Initially  Bent  Beam  on  an  Elastic  Foundation 

The  initial  gap  between  a beam  and  elastic  foundation 
is  given  by  the  formula  G [ l - x_  ],  where  x is  the  dis- 
tance from  the  center  of  the  bea?n,  and  G is  a constant 
greater  than  zero.  Matrices  B,  A,  b,  and  C are  the  same 
as  for  the  problem  in  the  preceding  section,  except  that 
the  vector  a is  obtained  from  the  expression  given  above. 

The  quadratic  programming  problem  was  solved  for  25 
potential  contact  points  on  the  half  beam,  with  an  inter- 
val size  of  0.06  in.  The  potential  contact  length,  as 
shown  in  Figure  10(a)  is  2.88  in.  With  G = 0.0001  the  re- 
sulting contact  length  was  1.38  in.,  for  loads  1000  lbs. 
and  2000  lbs.  These  results  are  shown  in  Figure  10(b). 

The  quadratic  programming  code  ZORILLA  £l0]  was  again  used 
to  solve  this  problem. 
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TABLE  2 


CONTACT  STRESS  DISTRIBUTION  AND  FINAL  GAP 
FOR  BEAM  ON  ELASTIC  FOUNDATION 
t = 1000  lbs. , U = 0.01 


POINT 

NO. 

CONTACT 
STRESS 
lb./sq. in. 

TOLERANCE 
LIMIT 
(in. ) 

FINAL  GAP  SIZE 
FOR  MODIFIED 
CONTOURS  (in.) 

1 

562.886 

0. 1000000E-01 

0 . 9999964E-02 

2 

562.888 

0 . 9982634E-02 

0.9910598E-02 

3 

562.889 

0 . 9930551E-02 

0.9614620E-02 

4 

562.891 

0 . 9843744E-02 

0.9165816E-02 

5 

562.891 

0.9722218E-02 

0 . 8595929E-02 

6 

462.892 

0 . 9565968E-02 

0 . 7923093E-02 

7 

562.892 

0 . 9374995E-02 

0.7147454E-02 

8 

562.894 

0.9149302E-02 

0.6315298E-02 

9 

562.893 

0 . 8888885E-02 

0 . 5400788E-02 

10 

562.893 

0 . 8593746E-02 

0 . 4458770E-02 

11 

562.895 

0.82638863-02 

0 . 3493937E-02 

12 

562.888 

0 . 7899299E-02 

0 . 2511400E-02 

13 

562.891 

0 . 7499997E-02 

0. 151101 9E-02 

14 

562.885 

0 . 7065967E-02 

0. 5030558E-02 

15 

393.626 

0.6597217E-02 

0.0 

16 

59.206 

0 . 6093744E-02 

0.0 

17 

0.0 

0 . 5555548E-02 

0.0 

18 

0.0 

0.4982639E-02 

0.0 

19 

0.0 

0 . 4374996E-02 

0.0 

20 

0.0 

0 . 3732643E-02 

0.0 

21 

0.0 

0 . 3055555E-02 

0.0 

22 

0.0 

0 . 2343752E-02 

0.0 

23 

0.0 

0 . 1597222E-02 

0.0 

24 

0.0 

0.8159771E-03 

0.0 

25 

0.0 

0.0 

0.0 

Contact  Stress  (lbs./sq. 


Contact  Length  L (in.) 

(b) 

Figure  10.  Contact  Stress  vs.  Contact  Length  from  Quadratic 

Programming  Solution,  with  G = 0.0001. 
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One  now  has  the  contact  arc  to  be  employed  for  solv- 
ing the  linear  programming  problem: 

Minimize  — = -b  + l 

u n 

subject  to  constraints 

^ *"n  + 1 ^ 0 

-b  £ -b° 
b ^ b1 


A A A A A 

-BS  - Aq  - b = a 


A 

A 


T A 

s 


C 


(3.2.1) 


VA  a/ 

-BS  - Aq 


£ 


A 

S > 0 


a + b 


A linear  programming  tableau  is  formed  from  Equations 
3.2.1,  as  explained  in  the  previous  section.  The  linear 
programming  problem  was  solved  and  results  are  presented 
in  Figures  (11-16).  A comparison  of  the  results  for  two 
loads,  t = 1000  lbs.  and  t = 2000  lbs.,  is  given  in 
Table  3.  Values  of  U used  are  proportional  to  applied 
load  for  the  cases  t = 1000  lbs.  and  t = 2000  lbs.  The 
peak  stress  obtained  for  t = 1000  lbs.  is  exactly  half 
that  for  t = 2000  lbs.  and  the  contact  length  is  the 
same  for  t = 1000  lbs.  and  t = 2000  lbs.  Thus,  contact 
arc  length  and  peak  stress  are  dependent  on  both  value 
of  applied  load  and  U.  Table  4 gives  the  optimum  stress 


Contact  Stress  (lbs./sq.  in. 


9 


Contact  Length  (in.) 


Figure  11.  Contact  Stress  vs.  Contact  Length  with  U = .01  and  G = 0.0001 
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Potential  Contact  Point 


Figure  12.  Gap  Size  vs.  Potential  Contact  Point  with  U = 0.01  and  G = 0.0001. 
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Figure  14,  Gap  Size  vs.  Potential  Contact  Point 

with  G = 0.0001. 
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Figure  15. 


Contact  Stress  vs. 

and 


Contact  Length  with  Augmented  Cost  Function 
with  G = 0.0001. 
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Figure  16.  Gap  Size  vs.  Potential  Contact  Point  with  Augmented  Cost 

Function  and  with  G = 0.0001. 
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TABLE  3 


COMPARISON  OF  RESULTS  FOR  LOADS  t = 1000  lbs. 
AND  t = 2000  lbs.  FOR  AN  INITIALLY  BENT 
BEAM  ON  ELASTIC  FOUNDATION 


LOAD  SR. 

U 

PEAK 

COMPUTING 

NO. OF  TIME/ 

CONTACT 

t NO. 

STRESS 

TIME 

ITRNS  ITRN 

LENGTH 

lbs./sq.in.  (SEC.) 

(SEC. ) 

(in.  ) 

1 

0.005 

677.6 

21.90 

3 

7.3 

1.62 

2 

0.01 

561.7 

39.05 

5 

7.81 

1.86 

1000 

3 

0.015 

498.7 

49.32 

6 

8.22 

1.98 

!• 

4 

0.02 

457.1 

71.82 

8 

8.98 

2.22 

1 

0.01 

1356.0 

21.65 

3 

7.217 

1.62 

2 

0.02 

1124.05 

41.37 

5 

8.274 

1.86 

2000 

3 

0.03 

997.2 

49.17 

6 

8.195 

1.98 

4 

0.04 

914.6 

73.74 

8 

9.22 

2.22 
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TABLE  4 

CONTACT  STRESS  DISTRIBUTION  AND  FINAL  GAP  FOR 
INITIALLY  BENT  BEAM  ON  ELASTIC  FOUNDATION 
t = 1000  lbs. , U = 0.01 


POINT  CONTACT  TOLERANCE 


INITIAL 


FINAL 


NO 

STRESS  LIMIT 

lb./sq.  (in.) 

in. 

GAP 
(in.  ) 

GAP 
(in.  ) 

1 

561.712 

0 . 1000000E-01 

. 1000000E-03 

0 . 9999931E-02 

2 

561.714 

0 . 9982634E-02 

0 . 9982637E-04 

0 . 9929962E-02 

3 

561.717 

0.9930551E-02 

0 . 9930554E-04 

0 . 9638943E-02 

4 

561.719 

0 . 9843744E-02 

0.9843748E-04 

0.91 967 58E-02 

5 

561.721 

0 . 9722218E-02 

0 . 9722220E-04 

0 . 8627117E-02 

6 

561.721 

0 . 9565968E-02 

0 . 9565971E-04 

0.7991239E-02 

7 

561.720 

0. 9374995E-02 

0 . 9374999E-04 

0 . 7174265E-02 

8 

561-720 

0 . 9149302E-02 

0 . 9149304E-04 

0 . 6328776E-02 

9 

561.722 

0 . 8888885E-02 

0 . 8888886E-04 

0 . 5433228E-02 

10 

561.722 

0 . 8593746E-02 

0. 8593748E-04 

0 . 44 8 98 2 8 E- 02 

11 

561.722 

0 . 8263886E-02 

0. 8263886E-04 

0 . 3522005E-02 

12 

561.722 

0. 7899299E-02 

0. 7899302E-04 

0 . 253201 6E-02 

13 

561.717 

0. 7499997E-02 

0 . 7499997E-04 

0.1539212E-02 

14 

561.719 

0 . 7065967E-02 

0 . 7065968E-04 

0.5400113E-02 

15 

406.454 

0.6597217E-02 

0.6597218E-04 

0.0 

16 

62.750 

0 . 6093744E-02 

0 . 6093745E-04 

0.0 

17 

0.0 

0 . 5555548E-02 

0. 5555548E-04 

0.0 

18 

0.0 

0 . 4982639E-02 

0 . 4982639E-04 

0.0 

19 

0.0 

0 . 4374996E-02 

0 . 4982639E-04 

0.0 

20 

0.0 

0 . 3732644E-02 

0 . 3732644E-04 

0.0 

21 

0.0 

0 . 3055555E-02 

0 . 3055555E-04 

0.0 

22 

0.0 

0 . 2343752E-02 

0 . 2343752E-04 

0.0 

23 

0.0 

0 . 1597222E-02 

0 . 1597222E-04 

0.0 

24 

0.0 

0.8159771E-03 

0.8159771E-05 

0.0 

25 

0.0 

0.0 

0.0 

0.0 
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distribution,  tolerance  limit,  initial  gap,  and  final 
gap  at  the  potential  contact  points.  Peak  contact  stress 
with  U = 0.01  and  t = 1000  lbs.  was  found  to  be  561.712 
lbs./sq.  in.  This  peak  stress  is  not  much  different  from 
the  one  obtained  in  the  case  of  beam  on  elastic  founda- 
tion with  no  initial  gap. 

The  quadratic  programming  problem  was  again  solved 
for  the  same  data,  but  with  an  initially  bent  beam  on  an 
elastic  foundation  with  G = 0.001  and  G = 0.005.  The  re- 
sults are  presented  in  Figures  (17-18).  It  is  observed 
that,  as  the  value  of  G is  increased,  peak  contact  stress 
decreases.  In  the  case  of  G = 0.005  contact  length  for 
1000  lbs.  and  2000  lbs.  loads  differs.  It  is  also  ob- 
served that  peak  contact  stress  depends  on  the  applied 
load . 

Once  the  contact  arc  is  known  from  the  contact  analy- 
sis problem,  the  linear  programming  problem  was  solved 
for  G = 0.001  and  G = 0.005,  with  two  different  loads  in 
each  case.  Results  are  presented  in  Figures  (19-22). 

Table  5 shows  the  comparison  of  results  with  U = 0.01, 
for  a beam  on  an  elastic  foundation  with  no  initial  gap 
and  with  an  initial  gap.  Table  6 shows  the  comparison  of 
results  for  t = 2000  lbs.  with  values  of  U ranging  from 
0.01  to  0.05.  It  is  noted  that  the  peak  contact  stress 
decreases  if  the  value  of  U and  G are  increased. 


Contact  Stress  (lbs./sq.  in. 


Contact  Length  (in.) 


Figure  17.  Contact  Stress  vs.  Contact-Length  from  Quadratic 

Programming  Solution  with  G = 0.001. 


fO 


Contact  Length  (in.) 


Figure  18.  Contact  Stress  vs.  Contact  Length  from  Quadratic 
Programming  Solution  with  G = 0.005. 


u> 


Contact  Length  (in.) 

Figure  19.  Contact  Stress  vs.  Contact  Length  with  G = 0.001. 


Potential  Contact  Point 

Figure  20.  Gap  Size  vs.  Potential  Contact  Point  with  G 


Contact  Stress  (lbs./sq.  in. 


Contact  Length  (in.) 


Figure  21.  Contact  Stress  vs.  Contact  Length  with  G = 0.005. 


cn 


V 


Potential  Contact  Point 


Figure  22.  Gap  Size  vs.  Potential  Contact  Point  with  G = 0.005. 


'j 
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TABLE  5 

COMPARISON  OF  RESULTS  FOR  BEAM  ON  ELASTIC  FOUNDATION 
WITH  NO  INITIAL  GAP  AND  WITH  DIFFERENT 
INITIAL  GAPS  FOR  U = 0.01 


LOAD 
(lbs. ) 

G 

CONTACT 
STRESS 
(lbs./ 
sq.  in.) 

CONTACT 
LENGTH 
(in. ) 

NUMBER 

OF 

ITRATIONS 

TIME/ITRN 
(sec. ) 

0 

562.9 

1.86 

5 

7.36 

1000 

.0001 

561.7 

1.86 

5 

7.81 

.001 

556.4 

1.86 

4 

8.6 

.005 

533.4 

1.98 

4 

9.375 

0 

1358.5 

1.62 

3 

7.25 

.0001 

1356.0 

1.62 

3 

7.22 

2000 

.001 

1349.5 

1.62 

2 

7.92 

.005 

1317.8 

1.62 

2 

8.5 

TABLE  6 


COMPARISON 

OF  RESULTS  FOR 

BEAM  ON 

ELASTIC 

FOUNDATION 

WITH  NO 

INITIAL  GAP  AND 

WITH  DIFFERENT 

INITIAL 

GAPS  FOR  t 

= 2000 

lbs . 

G 

ITEM 

QP 

LINEAR  PROGRAMMING 

SOLUTION 

SOLUTION 

U 

0.01 

0.02 

0.03 

0.04 

0.05 

0 

2242 

1358.5 

1125.6 

993.25 

914.8 

852.87 

.0001 

CONTACT 

2240 

1356 

1124 

997.2 

914.6 

852.6 

.001 

STRESS 

2220 

1349.5 

1118.4 

994.5 

911.1 

849.2 

.005 

2 

lbs. /in. 

2170 

1317.8 

1094 

973.3 

895.2 

833.9 

0 

1.38 

1.62 

1.86 

1.98 

2.22 

2.34 

.0001 

CONTACT 

1.38 

1.62 

1.86 

1.98 

2.22 

2.34 

.001 

LENGTH 

1.50 

1.62 

1.86 

1.98 

2.22 

2.34 

.005 

in. 

1.50 

1.62 

1.86 

2.10 

2.22 

2.34 
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Chapter  IV 
CONCLUSIONS 


The  technique  presented  here  to  solve  surface  con- 
tour design  problem  is  quite  simple,  and  has  given  con- 
sistently good  results.  By  adjusting  the  contours  of 
the  contacting  bodies,  the  stresses  developed  are  re- 
duced considerably.  The  problem  of  a beam  on  an  elastic 
foundation  demonstrates  some  of  these  facts.  The  peak 
stress  of  the  unmodified  structure  is  highest  at  the 
center  of  the  beam  and  decreases  at  the  end  points  of 
the  contact  region.  The  contact  design  problem  gives 
a relatively  constant  stress  distribution  on  the  con- 
tact region,  with  much  reduced  peak  contact  stresses. 

Tolerance  limits  affect  the  solution  as  is  shown 
in  Table  6.  In  this  problem,  the  load  and 
design  tolerance  limits  have  a measurable  effect  on  the 
solution  of  the  design  problem.  However,  as  shown  in 
Table  1,  the  solution  changes  proportionally  if  the 
load  and  tolerance  limit  are  changed  in  the  same 
proportion. 
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Suggestions 

This  design  method  can  be  extended  to  other  prob- 


lems  as 

follows : 

1. 

Asymmetrical  problems  for  beam  on  an  elastic 
foundation. 

2. 

Circular  inclusion  problems  for  elastic  bodies 
in  contact. 

3. 

Multibody  contact  problems. 
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APPENDIX  As 


FORMULATION  OF  MATRICES  B,  A,  a,  and  C 


!• 
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The  compatibility  condition  for  deformation  is 
given  by 

£ = u1  + u2  + d :>  0 (A-l) 

where  £ is  the  gap  vector  after  deformation  of  two 
. 1 2 

bodies,  u and  u denote  the  normal  displacement  vec- 
tors of  potential  contact  points  on  Bodies  1 and  2, 
and  d is  the  vector  of  initial  gap  between  the  two 
bodies. 

Body  2 of  Figure  23  is  fixed,  so  only  rigid 
body  displacement  for  Body  1 is  considered.  The  dis- 
placement of  points  on  Body  l is  due  to  rigid  body 
displacements  and  elastic  deformation.  Elastic  de- 
formation of  points  on  Body  l is  determined  relative 
to  zero  values  of  the  rigid  body  coordinates,  and 
total  displacement  is  found  by  superposition.  Gen- 
eralized displacement  coordinates  are  denoted  by  q. 

The  total  displacement  of  potential  contact 
points  on  Body  1 can  be  written  in  vector  form  as 

u1  = + Aq  (A- 2) 

e 

where  u1^  is  a reduced  vector  of  elastic  displacements 
of  potential  contact  points  on  Body  l.  The  vector  u1e 
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Body  2. 


S1  S2S3  S. 


2 2 2 
u i u 2 u 3 u2 


(b) 


Figure  23.  Two  Bodies  in  Contact. 
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is  determined  with  q = 0,  so  it  cannot  contain  dis- 
placement components  that  are  used  as  rigid  body  de- 
grees of  freedom.  The  projection  matrix  P is  given 
by  P = [0  | In-m-^'  w^ere  n t^ie  number  of  potential 
contact  points  and  m is  the  number  of  degrees  of  free- 
dom. The  matrix  A is  an  nxm  affine  transformation 
that  contributes  to  rigid  body  displacements. 

The  deformation  vector  uXe  is  further  decomposed 
as 

u1e  = F1PS  + V1e  (A-3) 

where  F1  is  the  flexibility  matrix  for  Body  1,  with 
q = 0,  and  VXe  is  the  vector  of  elastic  displacements 
of  potential  contact  points  due  to  externally  applied 
forces  tX  on  Body  1,  with  q = 0.  From  equations  (A-2) 
and  (A-3), 

U1  = PT[F1PS  + VXe]  + Aq  (A-4) 

For  Body  2, 

U2  = F2S  + V2e,  (A-5) 

2 . ... 

where  F is  an  nxn  flexibility  matrix  for  potential 

2 

contact  points  on  Body  2,  and  V is  displacement  of 
potential  contact  points  on  Body  2,  due  to  externally 
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applied  forces. 

Equation  (A-l)  can  now  be  written  as 
£ = (PTF1P  + F2)S  + Aq  + PTV1Q  + V2e  + d ^ 0,  (A-6) 


so 


T 1 2 

B = P F P + F 

A = A 


and 

T 1 2 

a = P VA  + V + d 
e e 

Equilibrium  Equations 

The  work  done  by  stresses  and  the  applied  forces 
on  Body  l,  due  to  rigid  body  displacements  only,  is 

T 

W = STAq  + t1  Hq  (A-7) 

where  H is  a matrix  that  gives  rigid  body  displacement 
at  the  points  of  application  of  t1 . 

Varying  the  rigid  body  displacements,  the  prin- 
ciple of  virtual  work  gives 

T 

SW  = STA5q  + t1  HSq  = 


0 


(A-8) 
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Thus, 


1 

(STA  + t1  H)5q  = 0 


(A-9) 


Since  all  components  of  q are  independent,  it  is 
necessary  that 


ATS  + HTt1  = 0 


(A-10) 


or 


T T 1 

AS  = -Ht 


Defining 


T 1 

C = -Ht 


(A-ll) 


Equation  (A-10)  can  be  written  as 


T 

AS  = C 


(A-12) 
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****************************************************** 
* * 

* MAIN  PROGRAM  * 

* BEAM  ON  AN  ELASTIC  FOUNOATION  WITH  NO  INITIAL  GAP.  * 

* * 
****************************************************** 


DESCRIPTION  OF  PARAMETERS. 

DEL  IS  THE  INTERVAL  SIZE  BETWEEN  THE  POTENTIAL  CONTACT 

POINTS ( P.C.P. ) 

AL  IS  HALF  THE  P.C. LENGTH. 

MNM  IS  THE  NO.  OF  DEGREES  OF  FREEDOM. 

NN  IS  THE  NO.  OF  P.C. PS. 

NC=NCI  IS  THE  NO.  OF  POINTS  IN  CONTACT. 

NNC  IS  THE  NO.  OF  POINTS  NOT  IN  CONTACT. 

P IS  THE  EXTERNALLY  APPLIED  LOAD. 

I TRN  IS  THE  LINEAR  PROGRAMMING  ITERATION  NUMBER. 

Ml  IS  THE  NO.  OF  CONSTRAINTS  PLUS  ONE. 

N 1 NO.  OF  DESIGN  VARIABLES  PLUS  NO. OF  GE.  TYPE  OF  CONSTRAINTS. 

M1=I+4*NC+MNM+NNC 

N1=2*MNM+1+2*NC+NGE 

DESCRIPTION  OF  MATRICES. 

A 1 ( NN* MNM)  IS  THE  AFFINE  TRANSFORMATION  MATRIX  FOR  RIGID  BODY 

DISPLACEMENT  CONTRIBUTION  FOR  BODY  1. 

T (NNf NN)  IS  THE  COMBINED  INFLUENCE  COEFFICIENT  MATRIX. 

F(NN,NN)  IS  EQUIVALENT  TO  T(NN,NN). 

FOLLOWING  MATRICES  ARE  GENRATED  FROM  THE  ABOVE  TWO  MATRICES. 


o> 

u> 


oooooooooooooooooooooooooonooooooo 


BHAT (NNyNN ) 
AHAT(NN»MNM) 
BTEL ( NN»NN ) 

ATEL(NN»MNM) 
B ( Ml yNl ) 


IS  THE  REDUCED  MATRIX  «T*  WITH  ROWS  AND  COLUMNS 
CORROSPONDING  TO  POINTS  NOT  IN  CONTACT  REMOVED. 

IS  THE  MATRIX  »AI'  WITH  CORROSPONDING  ROWS  FOR 
POINTS  NOT  IN  CONTACT  REMOVED. 

IS  THE  MATRIX  *T*  WITH  ROWS  CORROSPONDING  TO  POINTS 
IN  CONTACT  AND  COLUMNS  CORROSPONDING  TO  POINTS  NOT 
IN  CONTACT  REMOVED. 

IS  THE  MATRIX  »Al*  WITH  ROWS  CORROSPONDING  TO 
POINTS  IN  CONTACT  REMOVED. 

IS  THE  MATRIX  OF  CONSTRAINT  COEFFICIENTS. 


DESCRIPTION  OF  VECTORS. 


SA(NN) 

SSA(NN) 

SB(NN) 
SAHAT ( NN ) 

ASHAT ( NN) 

ASHAT(NN) 

SATEL(NN) 

SBTEL(NN) 

ASTEL(NN) 

RC( Ml ) 

NBP(Nl) 

Kll(NN) 

K (NN) 


IS  THE  VECTOR  OF  TOLRANCE  LIMIT  FOR  CONTOUR 
DESIGN  VARIABLES. 

IS  THE  VECTOR  OF  INITIAL  GAP  BETWEEN  THE  TWO  BODIES 
AT  THEIR  P.C.PS. 

IS  THE  VECTOR  OF  CONTOUR  MODIFICATIONS. 

IS  THE  VECTOR  * SA  * WITH  COMPONENTS  CORROSPONDING  TO 
POINTS  NOT  IN  CONTACT  REMOVED. 

IS  THE  VECTOR  * SSA ' WITH  COMPONENTS  CORROSPONDING  TO 
POINTS  NOT  IN  CONTACT  REMOVED. 

IS  THE  VECTOR  'SB*  WITH  COMPONENTS  CORROSPONDING  TO 
POINTS  NOT  IN  CONTACT  REMOVED. 

IS  THE  VECTOR  *SA»  WITH  COMPONENTS  CORROSPONDING  TO 
POINTS  NOT  IN  CONTACT. 

IS  THE  VECTOR  'SA*  WITH  COMPONENTS  CORROSPONDING  TO 
POINTS  NOT  IN  CONTACT. 

IS  THE  VECTOR*  SSA ' WITH  COMPONENTS  CORROSPONDING  TO 
POINTS  NOT  IN  CONTACT. 

IS  THE  VECTOR  OF  R.H.S.  VALUES  OF  CONSTRAINTS. 

COST  FUNCTION  ROW  VECTOR. 

IS  THE  VECTOR  OF  INDICED  FOR  POINTS  NOT  IN  CONTACT. 
IS  THE  VECTOR  OF  INDICES  OF  POINTS  IN  CONTACT. 


oooooooooooooooooooooooooooooooooo 


KIl(NN) 

IKNNI 

SHAT(NN) 
RBD(MNM) 
RBD1 ( 2*MNM ) 
7 ( MNM ) 

X ( I ) 

CST(NN) 

S IGN( Ml  ) 


IS  EQUIVALENT  TO  K(NN). 

IS  THE  VECTOR  HAVING  COMPONENTS  AS  I OR  0. 

1 FOR  POINTS  IN  CONTACT  AND  0 FOR  PTS.  NOT  IN  CONTACT. 
IS  THE  VECTOR  OF  CONTACT  FORCES. 

IS  THE  RIGID  BODY  DISPLACEMENT  VECTOR. 

IS  THE  VECTOR  HAVING  COMPONENTS  AS  TOTAL  EXTERNAL 
APPLIED  LOAD. 

IS  THE  DISTANCE  OF  I TH  P.C.P.  FROM  CENTER. 

IS  THE  CONTACT  STRESS  VECTOR. 

INDICATES  THE  TYPE  OF  CONSTRAINT  USED. 

-1  FOR  GE. TYPE  OF  CONSTRAINT. 

0 FOR  EQ.  TYPE  OF  CONSTRAINT. 

1 FOR  LE.TYPE  OF  CONSTRAINT. 


NOTE  THAT 


THE  DIMENSIONS  OF  THE  MATRICES  AND  VECTORS  SHOULD  BE  GIVEN  AS 
SHOWN  IN  THE  BRACKETS. ACTUAL  SIZE  OF  'HAT'  AND  'TELDA'  MATRICES 
AND  VECTORS  IS  GENERALLY  LESS. 

HAT  AND  TEL  ATTACHED  TO  THE  MATRICES  AND  VECTORS 
DENOTES  A AND  ~ RESPECTIVELY. 


INPUT  OATA  FOR  MAIN  PROGRAM  OF  LP  SUBROUTINE. 


1ST  SET  OF  CARDS. 
2ND  SET  OF  CARDS. 
3RD  SET  OF  CARDS. 
4TH  SET  OF  CARDS. 


READ  IN  MI,N1 

READ  IN  OR  GENERATE  MATRIX  *B* 

READ  IN  OR  GENERATE  VECTOR  *RQ* . 

READ  IN  SIGN  ACCORDING  TO  TYPE  OF  CONSTRAINT. 


IMPORTENT  NOTE. 


ooononoonnonooooooo 


IN  CASE  OF  BEAM  ON  ELASTIC  FOUNDATION  WITH  INITIAL  GAP 
A FEW  CHANGES  IN  THE  MAIN  PROGRAM  ARE  NECESSARY. 

VECTOR  SSA  IS  CHANGED. 

SSA  IN  THE  EXAMPLE  SOLVED  IN  THIS  WORK  IS  TAKEN  AS  GIVEN  BELOW. 

SSA< I) =0.0001* (l.-( (X< I))**2)/(AL**2)» 

FEW  CHANGES  IN  MATRIX  ‘B*  ARE  AS  FOLLOWS. 

B( I+3*NC,K1)=-AHAT< I,K1) 

320  B ( I + 3<cNC»  Kl+MNM ) = AHAT  ( I »K  1 ) 

B(I+3*NC,MM+I )=-l. 

330  B(I+3*NC, MM+NC+K1 ) =-BHAT ( I , K1 ) 


INTEGER  SIGNtGP*  PG 
REAL  NBP 
C 

COMMON  B( 105, 105) ,RQ( 105 ) ,NBP ( 105 ) , SIGN ( 105) * KK( 105 ) , I TRN 
COMMON  RRQ( 105  » ,KN,GP,PG 
C 

DIMENSION  SHAT (25), SA(25),ASHAT(25),SB(25),K(25), 11(25) 

DIMENSION  KI1(25) ,BHAT( 25,25 ),SAHAT( 25), AHAT( 25,2) , SSA (25) 
DIMENSION  BTEL(25,25» , ATEL( 25 , 2 ) , SATEL ( 25 ) , SBTEL ( 25 ) ,ASTEL(25» 
DIMENSION  Kll ( 25 ) , P( 105 ) , RBD( 2 ) ,RBDl ( 4 ) 

DIMENSION  T(25,25»,A1(25, 2),Z( 2),SBHAT( 25 ) , F( 25, 25 ) , X ( 25 ) ,CST( 25 ) 
C 

READ(5,3)  PP,U,AL,DEL,MNM,NN,NC,NNC,NGE 
C 

3 FORMATUFIO.A,  515) 


o o o o 


z<n=pp 

Z ( 2 ) =0. 0 
DO  5 J=1,MNM 
C 

5 READ(5,2)(A1(I,J),I=1,NN) 

C 

DC  20  1=1, NN 
20  S SA ( I ) =0 . 

DC  1 1=1, NN 
X ( I ) = ( I -1 ) *DEL 

1 SA( I )=U*( 1 .-( ( t X { I ) >**2)/(AL**2) )) 

DC  4 1=1, NN 

SB ( I » =0 . 

C 

4 READ(5,2)(F(I,J),J=1,NN) 

C 

2 FORMATt 5E15.7) 

DC  22  1=1, NN 
DC  22  J=1 , NN 

22  T (I,J)=F(I,J) 

I TRN= 1 

K(NC)  EQUIVALENT  TO  KIKNCI  IS  THE  INDEX  OF  POINTS  IN  CONTACT. 


READ( 5,40) (K( I ) , 1=1, NC) 
C 

40  FORMAT (401 2) 

DO  25  1=1, NC 
25  KI1 ( I )=K( I ) 

725  CONTINUE 
JJ=0 

DC  35  1=1, NC 
35  K ( I ) =KI  1 ( I ) 


<! 


d 


K0=0 
NM=0 
K2=0 
N 12=0 

DC  42  1=1, NN 
4?  11(11=0 

DC  60  1=1, NN 
DC  70  NN1=1,NC 
L=K ( NN1 ) 

I F ( I-L ) 69, 100,69 
100  KC=K0+1 
I I ( L ) =1 
C 

C TG  FORM  MATRIX  B-» HAT • , A- 1 HAT • , SMALL-A- * HAT* , SMALL-B- «HAT • , 

C AND  VECTOR  OF  INITIAL  GAP  FOR  POINTS  IN  CONTACT. 

C 

DO  220  KUI=1,MNM 
220  AHAT  ( KO » KU I ) = A 1 ( L , KU I I 
S AHAT (KO ) =SA ( L ) 

SBHAT ( KO) =SB ( L ) 

ASHAT(KO)=SSA(L) 

DC  71  J= 1 , NN 
DC  71  KM= I , NC 
M3=K{ KM ) 

I F ( J-M3)71,200,71 
200  NM=NM+ 1 

BHAT(KO,NM)=T( I, J) 

71  CONTINUE 
NM=0 

GC  TO  70 
69  K2=K2+ 1 

IF(NC~K2)70,80,70 
80  JJ=JJ+1 


CTi 

09 


oooo  oonn 


TO  FORM  MATRIX  B-TELDAI BTEL1 , A-TELDA ( ATEL I , SMALL-A-TELDA 
(ATEL),SMALL-B-TELDA(BTEL>,ANO  VECTOR  OF  INITIAL  GAP  FOR 
POINTS  NOT  IN  CONTACT. 

DO  222  MZ  = I , MNM 
222  ATEU  JJ,MZ)=A1(  I , MZ  I 
$ATEL< JJ>=SA<  I ) 

. SBTEL( J J)=SB(  I ) 

ASTELI JJ)=SSA< I ) 

K 1 1 ( JJ ) =1 
DC  72  J = 1 , NN 
DC  72  KM3= 1 , NC 
M2=K( KM3 ) 

I F ( J-M2)72,250,72 
250  N 12=N12+ 1 

BTELI JJ,N12)=T( I,J> 

72  CCNTINUE 
N 12  = 0 

70  CCNTINUE 
K2  = 0 

60  CCNTINUE 
101  F CRM AT (/5E15.7/J 

N1=2*MNM+1+2*NC+NGE 
N=2*MNM+1+2*NC 
M1=1+4*NC+MNM+NNC 

FCRMULATION  OF  *B'  MATRIX. 

FCRMUIATI ON  OF  VECTOR  • RQ • . 

DO  300  1=1, Ml 
RC( I )=0. 

DC  300  J=1,N1 
300  B(I,J)=0. 

C 


C COST  FUNCTION  ROW  FOR  MAXIMIZATION  PROBLEM(NBP) . 

c 

DO  310  1=1, N1 
310  NBP(NGE+I)=0. 

MM=2*MNM+ 1 
NBP(NGE+MM)=-1 .0 
DC  350  1=1, NC 
B(I,MM)=-i.O 
B ( I , MM+ NC+ 1 1=1.0 
B ( I +NC , MM+ I )=1. 

B(I+2*NC,MM+I)=-1. 

DO  320  Kl= 1 , MNM 
B(I+3*NC,K1>=AHAT< I,Ki) 

320  B(I+3*NC»K1+MNM) =-AHAT ( I , K1 ) 

B(I+3*NC,MM+I)=1. 

DO  330  Kl= 1 » NC 

330  B ( I +3*NC , MM+NC+K1 I «BHAT ( I ,K 1 ) 

DC  340  Ki= 1 , MNM 

340  8<4*NC+Kl,MM+NC+I)=AHAT< I ,K1) 

350  CONTINUE 

DO  360  KJ= 1 , NNC 
DC  355  Kl= 1 , MNM 
B(KJ+4*NC+MNM,K1 ) =-AT EL < K J , K 1 ) 

355  B(KJ+4*NC+MNM,K1+MNM) =ATEL(KJ,K1> 

DC  358  KU= 1 , NC 

358  B ( KJ+4*NC+MNM, MM+NC+KU I =-BTEL ( KJ ,KU ) 

360  CONTINUE 

DO  370  J=I , NC 

RC( J ) =0 .0 

RQ( J+NC )=SAHAT ( J) 

RC( J+2*NC)=SAHAT( J) 

370  RC< J+3*NC»=ASHAT( J) 

DO  375  11=1, MNM 
375  RC( 4*NC+1 1 ) =Z ( 1 1 ) 


non  n noon  nooooo 


V 


00  380  IJ  =1 » NNC 

380  RC(4*NC+MNM+IJ)=SBTEL< IJ»+ASTEL< IJ» 

SIGN  INDICATES  THE  TYPE  OF  CONSTRAINTS. 

ONE  FOR  LESS  THAN  EQUAL  TO  TYPE  OF  CONSTRAINT. 

ZERO  FOR  EQUALITY  CONSTRAINT. 

MINUS  ONE  FOR  GREATER  THAN  EQUAL  TYPE  OF  CONSTRAINT. 

S IGN (Ml >=0 
MK=3*NC 
DC  399  1=1, MK 
399  S IGN( I ) =1 
KM1 =MK+ 1 
KM2=MK+NC+MNM 
DC  420  I =KM1 » KM2 
420  S IGN ( I ) =0 
KM3=KM2+1 
KM4=KM2+NNC 
DC  430  I=KM3»KM4 
430  S IGN ( I » = 1 

ALL  THE  COEFFICIENT  IN  THE  LP  TABLEAU  START  LEAVING  NGE- 
CCLUMNS  BLANK. 

DC  600  1=1, Ml 
DC  600  J=1,N 
600  B( I ,NGE+J)=B( I , J) 

GP=0 

PG=0 


CALL  LINP  (M1,N,NGE) 

CHECK  IF  SOLUTION  IS  UNBOUNDED  OR  INFEASIBLE  THEN  TERMINATE. 


'O 

H* 


IF(GP.EO.l)  GO  TO  1000 
IF(PG.EQ.l)  GO  TO  1000 
00  308  J=1 ,KN 
C 

C REARRANGING  THE  INDICIES  OF  THE  VARIABPES  IN  A SERIAL  OROER. 

C 

DC  308  1=2, KN 
P(J)=KK( J) 

IF< J-I 1307,308,308 

307  IF(P(J)-KK( I ) >308,308,299 
299  KK ( J ) =KK ( I ) 

KK( I ) =P ( J ) 

308  CONTINUE 

DO  298  1=1, KN 
298  RC(KK( I ) )*RRQ( KK (111 
C 

C RBD  IS  RIGID  BODY  DISPLACEMENT. 

C 

I M=2*MNM 
DO  290  1 = 1, IM 
290  RBDI ( I ) =0. 

C 

C TO  FIND  THE  RIGID  BODY  DISPLACEMENTS  ( SMALL-*Q* ) 

C 

DC  205  1 = 1, IM 
I F ( KK ( I )-IM) 199, 199,205 
199  KX= I 

I F(  KK  ( I ) . L E . 2 ) GO  TO  202 
GO  TO  275 

202  I F ( KK ( I ) • E Q*  2 ) GO  TO  280 
GO  TO  285 

280  RBDl ( 2 ) =-RQ( KK ( I ) ) 

GO  TO  205 

285  R8D1( 1 )=RQ(KK( I ) ) 

-j 

to 


■ 
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GC  TO  205 

275  IFUM.EQ.2)  GO  TO  205 

282  IF(KK( I I.EQ.3I  GO  TO  283 
GC  TO  284 

283  RBD1(3>=RQ(KK< I » ) 

GC  TO  205 

284  I F( KM  I ) .EQ.4 ) GO  TO  279 
GO  TO  205 

279  RBD1 (4  I =-RQ( KK ( I ) ) 

205  CONTINUE 
N21=0 

DO  281  J = 1 » MNM 

RBD( J)=R8D1( J+N21J+RBDK J+1+N21) 

N21=N21+1 

281  WRITE(6,214)  N21 , RBO( J ) 

214  FORMAT < /2X, *RBD< • , 1 1 , • I = • , E15.7/ ) 

KX=KX+ 1 

FX2=RQ(KK(KXJ ) 

TO  GET  THE  CONTOUR  DESIGN  VARIABLES. 

CCNTOUR  DESIGN  VARIABLES  START  WITH  INDEX  OF  RQ  AS  2*MNM+2. 

J 1=NC+MM 
DC  210  1=1, NC 
J2=KX+1 
JC= I +MM 

IF(KK( J2I-J1 1206,206, 221 

206  CONTINUE 

IFIKM J2). EQ.JOI  GO  TO  203 
. GO  TO  221 
203  KX=KX+ 1 

SB(K( I ) )=RQ(KK( J2) ) 

GO  TO  210 

221  S8(K(in=0.0 


-j 

co 
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210  CONTINUE 

WRITE ( 6 » 226  > FX2 

226  FCRMAT(/2X»'C0ST  FN=',E15.7/» 

NCI=NC 
WRITE! 6,211) 

211  FORMAT </ IX, 'SR.NO. ',3X, 'POINT* ,9X, ' SHAT* , 15X, 'EPHAT' , 15X, 'SBHAT* ) 

CONTACT  STRESS  VECTOR  FORMATION. 

DO  224  1=1, NC 
JK=KX+l 
K 12=MM+NC+ I 

IF(KK( JK» .EQ.K12I  GO  TO  277 
GO  TO  278 

277  KX=KX+ 1 

SHAT ( I ) =RQ( K12 ) 

GC  TO  224 

278  SHAT (11=0.0 

224  CONTINUE 

COMPUTE  EPSLON  'HAT*. 

OC  225  1=1, NC 
EPHAT=SB(K( I ) )+ASHAT(  I ) 

DO  208  K8= i , MNM 

208  EPHAT=EPHAT+AHAT( I , K8 I *RBD( K8 > 

00  209  J= 1 , NC 

209  EPHAT=EPHAT  + BHAT( I , J ) *SHAT(  J ) 

225  WRITE(6,223)  I , K ( I ) , SHAT< II , EPHAT, SB( M I )) 

223  FORMAT (/5X, 12, 3X, 1 2, 3X , El 5 . 7 , 5X, E15.7, 5X, E15. 7) 

COMPUTE  EPSLON  'TELDA'. 

CHECK  WHICH  OF  THE  EPSLON  TELDA  IS  ZERO. 


o o o o o o 


WRI TE ( 6 ,235 ) 

235  FORMAT ( /5X , • PO INT*  * 8X » * EPSLON  TELDA*/) 

DC  230  1=1, NNC 

EPTL=S8TEL( I ) +ASTEL ( I ) 

DC  236  K7= 1 , MNM 

236  EPTL=EPTL+ATEL( I ,K7)*RBD(K7> 

DC  240J=1, NC 

240  EPTL=EPTL+BTEL< I , J)*SHAT( J) 

WRI T E ( 6 , 23 1 ) K 1 1 ( I ) , EPTL 
IF(EPTL.LE.lE-04)  GO  TO  232 
GC  TO  230 
232  II(Kll(in=l 

230  CCNTINUE 

231  FCRMAT(/5X, I5,5X,E15.7> 

S-HAT  IS  THE  VECTOR  OF  CONTACT  FORCES. 

CHECK  WHICH  OF  THE  S-HAT  ARE  ZERO. 

POINTS  WHERE  S-HAT  IS  ZERO  ,IS  LIFTED. 

KNC  ARE  THE  NO.  OF  PTS.  LIFTED. 

KNC  = 0 

251  FORMAT ( /2X, 'KNC-NO. OF  PTS.  WHICH  CAN  BE  LIFTED=',I3) 
DO  259  1=1, NC 
I F( SHAT ( I I 1259,249,259 
249  KNC=KNC+1 
I I ( K( I ) ) =0 
259  CCNTINUE 

WRITE(6,251) KNC 
NC=0 

DO  270  1=1, NN 
I F ( 1 1 ( I ) . EQ. 0 ) GO  TO  270 
NC=NC+ 1 
K I 1 (NC ) =1 
270  CONTINUE 


m 


non 
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NNC=NN-NC 

WRITE(6,255)  NC,NNC 

255  FCRMAT(/2X,'NC=',I3,2X,  *NNC=' »I3/J 

COMPARE  THE  CONTOUR  BETWEEN  THE  TWO  CONSECTIVE  ITRNS. 

IF(NCl-NC)  725,800,725 
800  CONTINUE 
J IK  = 0 

DO  825  1=1, NC 

IF(K( IJ-KI 1( I» ) 728,825,728 
728  J IK  = J IK  + 1 
825  K C I J =KI 1 ( I ) 

IF(JIK.EQ.O)  GO  TO  700 
GO  TO  725 
700  WRI TE ( 6 , 860 ) 

860  FCRMATI/5X, 'FINAL  RESULTS'/) 

WRITE(6,875)  FX2 

875  FCRMATI/1X, 'COST  FN=',Fl2.6/> 

WRIT  E ( 6 ,405 ) 

405  FCRMATI/20X, • INITIAL  GAP'/) 

WRITE (6, 2) ( S S A C I ) , 1 = 1, NN) 

WRITE(6,403) 

403  FORMAT(/20X, 'TOLERANCE  LIMIT'/) 

WRITE(6,2) (SA( I ) , 1=1, NN) 

WRITE(6,880) 

880  FCRMATI/5X, 'POINT', 7X, 'FINAL  GAP  ',7X'C0NTACT  FORCE', 6X 
1 • CONT ACT  STRESS' ) 

OC  900  1=1, NC 
CST ( I ) =SHAT ( I ) /DEL 

900  WRI TE ( 6 ,890 ) K(I),SB(K(I)), SHAT ( I ) ,CST( I ) 

890  FORMAT (/5X,15,5X,E15»7,5X,E15.7,5X,E15.7) 

1000  CALL  EXIT 
END 


•vj 

CTi 


noon 


c ********************************** 

c * * 

C * SUBROUTINE  LINEAR  PROGRAMMING.  * 

C * * 

C ** ******************************** 

C 

SUBROUTINE  LINP  (M1,N,NGE) 

C 

INTEGER  RNM1 ♦ RNM2  tCLNMI , CLNM2 » BLNK  ,LPFTN  1 

I IBN1( 105) » IBN2( 105),NBN1( 105) ,NBN2( 105) 

C 

INTEGER  SIGN,GP,PG 

REAL  PIVOT,LST,XNBP»FN,CJBAR,X,VALUE,BP(105),PI(105) ,XPI(105) 

REAL  NBP 
C 

COMMON  B( 105, 105),RQ( 105 ) ,NBP ( 105 ) , SIGN ( 105 ) , KK ( 105 ) , I TRN 
COMMON  RRQ (105),KN,GP,PG 

C LPFTN  5 

OATA  BLNK/4H  / 

DATA  NM1 , NM2/ • CCCC ' , 1 AAAA ' / 

LPFTN  8 
LPFTN  16 

INPUT  PROGRAM  LPFTN  17 

LPFTN  18 

N I = 5 
NO  = 6 
I N=1 
M=M1-1 

DO  101  1=1, M 
IF(SIGN(I ) ) 108, 107, 106 
106  BP( I ) =0 . 

GO  TO  101 
108  BP( I ) =-l .0 
B ( I * IN) =-l .0 


•vl 

*■0 


o o o non 


NBNl ( IN)=NM2 

NBN2(IN)=IN 

NBP(IN»=0. 

I N= I N+ 1 
GC  TO  101 
107  BP ( I ) =-2.0 
101  CONTINUE 

DC  102  J = 1,N 
NBNl (J+NGE )=NM1 
102  NEN2( J+NGE )= J+NGE 
N=N+NGE 
DC  10  1 = 1, M 

I F( BP( I ) + l .0)  19,11,12 

11  I BN  1 ( I ) =BLNK 
I BN2 ( I ) =BLNK 
GC  TO  10 

19  B P ( I)=-l.O 
GC  TO  11 

12  CONTINUE 

I BN1 ( I ) =NM2 
I BN2 ( I ) = I 
10  CONTINUE 

LPFTN164 


ACCUMULATE  COUNT  OF  INFEASIBILITIES  LPFTN165 

LPFTN166 

NINF  =0  LPFTN167 

DO  6000  1=1, M LPFTN168 

I F ( BP ( I ) >6001,6000,6000  LPFTN169 

6001  NINF  = NINF+1  LPFTN170 

6000  CONTINUE 

LPFTN172 

GENERATE  INDICATORS  FOR  MINIMIZATION  OF  INFEASIBILITY  LPFTN173 

LPFTN174 
LPFTN175 


DO  6101  J= 1 , N 


-J 

CD 


ooo  nnnui  on 


XPI(J)  =0. 

DC  6101  1*1, M 
I F( BP( I I) 6102 ,6101, 61 01 
6102  XPI(J)  = XPI(J)-8(I,J) 

6101  CCNTINUE 

DC  6002  I = 1 » M 
6002  B P ( I ) = 0. 

WRITE( 6*401) 

401  FORMAT ( /2X , ' ************************* / ) 
WRI TE ( 6 1 400 ) I TRN 

400  FCRMAT(/5X,*MAIN  I TRN  N0.=»,I3/) 
WRITE(6,402) 

402  FORMAT ( /2Xt •************************•/) 
I TRN= I T RN+ 1 
I PHASE  = 1 

MAIN  ROUTINE 


LPFTN176 
LPFTN177 
LPFTN178 
LPFTN179 
LPFTN180 
LPFTN18 1 
LPF1N182 


LPFTN183 

LPFTN184 

LPFTN185 


9201  WRITEINO, 92021 

9202  FORMAT  («0  ITERATION  VAR  IN 
IT  = 0 

4325  CCNTINUE 

CALCULATE  SHADOW  PRICES 

DC  194  J= 1 » N 
PKJI  = -NBP(J) 

DC  194  1*1, M 

194  PI(J)  = PI(J)  + BPm*B(I,J) 

SELECT  BEST  NONBASIS  VECTOR 

9101  LST  = -.0000001 
KCOL  = 0 


VAR  OUT  OBJ  FN',/»  LPFTN188 

LPFTN189 

LPFTN190 

LPFTN191 

LPFTN192 

LPFTN193 

LPFTN194 

LPFTN195 

LPFTN196 

LPFTN197 

LPFTN198 

LPFTN199 

LPFTN200 

LPFTN202 


VO 


o o o n o o 


GO  TO  (751,5521, IPHASE 
751  IF(NINF)54321, 54321,55 2 
552  CONTINUE 

DG  9102  J=  l » N 

IGNORE  ARTIFICIAL  VARIABLES 

I F(NBN1 (J)-BLNK+NBN2( JJ-BLNK) 651,9102,651 
651  CONTINUE 

GC  TO  ( 6003,6004), IPHASE 

6003  I F( XPI ( J )-LST) 6005,6006,6006 

6005  KCOL= J 

LST  = XPIIJI 
GC  TO  9102 

6004  CONTINUE 

I F( PI (J)-LST) 9103,9102,9102 

9103  KCOL  = J 
LST  = PI ( J ) 

6006  CONTINUE 
9102  CONTINUE 

IF  (KC0D54321, 54321, 9104 

DETERMINE  KEYROW 

9104  KROW  = 0 
CJBAR  = LST 
LST  = 1.0E20 
DC  9105  1=1, M 

I F( B( I, KCOL) )9105, 9105, 9106 

9106  RATIO  = R0( I )/B( I, KCOL) 

IF  (RATI0-LSTI9107, 9105, 9105 

9107  LST  = RATIO 
KROW= I 

9105  CONTINUE 


LPFTN203 
LPFTN204 
LPFTN205 
LPFTN206 
LPFTN207 
LPFTN208 
LPFTN209 
LPFTN2 10 
LPFTN211 
LPFTN212 
LPFTN213 
LPFTN214 
LPFTN215 

LPFTN217 
LPFTN218 
LPFTN219 
LPFTN220 
LPFTN22 1 
LPFTN222 
LPFTN223 
LPFTN224 
LPFTN225 
LPFTN226 
LPFTN227 
LPFTN228 
LPFTN229 
LPFTN230 

LPFTN232 

LPFTN233 

LPFTN234 

LPFTN235 

LPFTN236 


— • 

• 

I F( KROW >9112,9112,9114 

LPFTN237 

9112 

WRITE (NO, 9 113)  NBN 1 ( KCOL) ,NBN2 ( KCOL ) 

9113 

FORMAT ( * VARIABLE  *,A2,I3,*  UNBOUNDED  *1 

GP=GP+1 

GO  TO  54323 

LPFTN240 

9114 

CONTINUE 

LPFTN241 

C 

LPFTN242 

C 

TRANSFORM 

LPFTN243 

C 

LPFTN244 

C 

DIVIDE  BY  PIVOT 

LPFTN245 

PIVOT  = B(KROW,KCOL) 

DC  9108  J=1,N 

LPFTN247 

9108 

B (KROW, J ) = B(KROW, JJ/PIVOT 

LPFTN248 

RC(KROW)  = RQ( KROW ) /P I VOT 

LPFTN249 

DO  9109  1=1, M 

LPFTN250 

I F( I -KROW  >9110,9109,9110 

LPFTN251 

9110 

RC(I)  = RQ( I ) - RQ(KROW)*B( I, KCOL) 

LPFTN252 

DC  4444  J = 1,  N 

I F ( J-KCOL) 91 11,4444, 91 11 

LPFTN254 

9111 

B ( I , J ) = B( I , J ) - B ( KROW, J ) *B ( I , KCOL ) 

LPFTN255 

4444 

CONTINUE 

9109 

CONTINUE 

LPFTN256 

DO  9300  1=1, M 

LPFTN257 

9300 

B ( I , KCOL ) = -B(I»KCOL)/PIVOT 

LPFTN258 

B ( KROW , KCOL ) = l.O/PIVOT 

LPFTN259 

C 

LPFTN260 

c 

INTERCHANGE  BASIS  AND  NONBASIS  VARIABLES 

LPFTN2 

c 

LPFTN262 

RNM1  = NBNKKCOL) 

LPFTN263 

RNM2  = NBN2( KCOL ) 

LPFTN264 

NBNKKCOL)  = I BN1  (KROW ) 

LPFTN265 

NBN2(KC0L)  = IBN2(KR0W) 

LPFTN266 

IBNl(KROW)  = RNMI 

LPFTN267 

I BN2 ( KROW ) = RNM2 

LPFTN268 

CD 

t— * 

o o o 


6200 

6201 

9301 

7000 

7003 

7001 

9302 

7004 
C 

6112 

6111 


LST  = NBP(KCOL) 

NBP(KCOL)  = BP(KROW) 

BP(KROW)  = LST 
IT  = IT  ♦ 1 

I F(  NBN1  (KC0L)-BLNK+NBN2{  KCOO-BLNK)  6201, 6200, 6201 

NINF  = NINF-i 

CONTINUE 


LPFTN269 

LPFTN270 

LPFTN271 

LPFTN272 

LPFTN273 

LPFTN274 

LPFTN275 


COMPUTE  OBJECTIVE  FUNCTION 

FN  = 0. 

DC  9301  1=1, M 

FN  = FN  ♦ BP( I ) *RQ( I ) 

GC  TO  (7000, 7001), IPHASE 
SAVE  = PKKCOL) 

DO  7003  J= 1 , N 

PIIJ)  = PI ( J ) - SAVE*B(KROW,J) 
XPI(J)  = XPI(J)  - CJBAR*B(KROW, J) 
CONTINUE 

PKKCOL)  = -SAVE/PIVOT 
XPI(KCOL)  = -CJBAR/P I VOT 
GC  TO  7004 
CONTINUE 
DC  9302  J= 1 , N 

PIIJI  = PI(J>  - CJBAR*B(KROW,J) 
PKKCOL)  = “CJ  BAR/PIVOT 
CONTINUE 

CHECK  FOR  ESSENTIAL  ZERO 
DC  6111  1 = 1, M 
DO  6111  J= 1, N 
X=B( I, J > 

I F ( ABS { X )- .0000001 )6112, 6112, 6111 
B(  I , J)  = 0. 

CONTINUE 


LPFTN277 
LPFTN278 
LPFTN279 
LPFTN280 
LPFTN28 1 
LPFTN282 
LPFTN283 
LPFTN284 
LPFTN285 
LPFTN286 
LPFTN287 
LPFTN288 
LPFTN289 
LPFTN290 

LPFTN292 

LPFTN293 

LPFTN294 

LPFTN295 

LPFTN296 

LPFTN297 

LPFTN298 

LPFTN299 

LPFTN300 

LPFTN301 

LPFTN302 


c 

c 

C LOG  ITERATION 

C 

WRITE (NO, 9120 1 IT, I BN1 < KROW ) , I BN2 (KROW) , NBN1 { KCOL) , NBN2 (KCOL ) ,FN 
9120  FORMAT ( 19, 7X, A2, 1 3, 8X , A2, I3,3X,F13.3» 

GO  TO  9101 
C 
C 

54321  CONTINUE 

I F( I PHASE- 1 18000,8000, 54322 
8000  I PHASE  = 2 

I F( NINF 18003,8003,8004 

8004  WRITEt  NO, 8005 ) 

8005  FORMAT ( ' 0 SOLUTION  INFEAS IBLE • , / I 
PG=  PG+1 

GO  TO  54323 
8003  CONTINUE 

WRITEt  NO, 8002 ) 

8002  FORMAT ( * 0 SOLUTION  FEASIBLE  ',/) 

GO  TO  54325 

54322  CONTINUE 
C 

C GUTPUT  ROUTINE 

C 

WRITE(N0,301»  IT,  FN 

301  FORMAT ( • 1 * , * ITERATION* , 15, • OBJ  FN  »,F15.3/) 

WRITEt  NO, 302 ) 

302  FORMATt 3X, 'BASIS  VAR ', 17X ,* AMOUNT ' ) 

KN=0 

DO  3033  1=1, M 
C 

C COST  RANGING 
C 


LPFTN303 
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LPFTN309 
LPFTN310 
LPFTN31 1 
LPFTN312 
LPFTN313 
LPFTN314 
LPFTN315 

LPFTN317 


LPFTN319 

LPFTN321 

LPFTN323 

LPFTN324 

LPFTN325 

LPFTN326 

LPFTN328 


LPFTN332 

LPFTN333 

LPFTN334 

LPFTN335 
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VALUE  = 1.0E20 
L ST  = 1.0E20 
DC  12300  J = 1 » N 

I F( NBNl ( J ) -BLNK+NBN2 ( J ) -BLNK ) 12305,12300,12305 
12305  CCNTINUE 

I F ( B ( I , J ) ) 12301,12300,12302 

12302  X=PI ( J)/B( I,J) 

I F(X-LST) 12303, 12300, 12300 

12303  LST=X 

GC  TO  12300 
12301  X=-PI(JJ/B(I,J) 

IF(X-VALUE)12304, 12300,12300 

12304  VALUE  = X 
12300  CCNTINUE 

LST  = BP( I>  - LST 
VALUE  = BP( I)  + VALUE 
C 

C VARIABLES  WITH  NAMES  AS  'CC»  ARE  SEPERATED  - 

C THESE  VARIABLES  ARE  THE  DESIGN  VARIABLES  OF  LP  PROBLEM. 

C 

C 

C VARIABLES  WITH  NAMES  ' AA'  ARE  ELIMINATED  IN  THE  NEXT  STEP. 

C THESE  ARE  THE  SLACK  VARIABLES. 

C 

I Ft IBN1( I J-NMl) 3033, 306,3033 
306  KN=KN+1 

KK(KN)=IBN2<  I ) 

RRQ(KMKN)  ) = RQ(  I ) 

WRITE (6, 304) IBN1( I ) , I BN2( I ) ,RQ( I ) 

3033  CCNTINUE 

304  FCRMAT(7X,A2, I3,7X,F16.6) 

54323  CCNTINUE 
RETURN 
END 


LPFTN336 

LPFTN337 

LPFTN338 

LPFTN339 

LPFTN340 

LPFTN341 

LPFTN342 

LPFTN343 

LPFTN344 

LPFTN345 

LPFTN346 

LPFTN347 

LPFTN348 

LPFTN349 

LPFTN350 

LPFTN351 


LPFTN  19 
LPFTN365 


CD 
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C *************** 
C * * 
C * INPUT  DATA.  * 
C * * 
C *************** 


p 

U 

AL 

DEL  MNM 

NN  NC  NNC 

NGE 

1000.  0 

.01 

1.44 

0.06 

1 

25  12 

13 

0 

Al- 

MATRIX. 

0.1000000E 

01 

0. 1000000E 

01 

0 . 1000000E 

01 

0. 1000000E 

01 

0. 1000000E 

01 

0. 1000000E 

01 

0. 1000000E 

01 

0. 1000000E 

01 

0.1000000E 

01 

0 • 1000000E 

01 

0 • 1000000E 

01 

0 • 1000000E 

01 

0 . 1000000E 

01 

0. 1000000E 

01 

0.1 OOOOOOE 

01 

0. 1000000E 

01 

0.  1000000E 

01 

0.1000000E 

01 

0. 1000000E 

01 

0. 1000000E 

01 

0.1000000E 

01 

0. 1000000E 

01 

0 . 1000000E 

01 

0. 100000CE 

01 

0. 1000000E 

01 

T-MATRIX. 


ROW  1 


0.5000000E-04 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 


O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 


O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 


O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 


O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 
O.OOOOOOOE  00 


ROW  2 

O.OOOOOOOE  00 
0. 1 42  3058E-06 
0.2947765E-06 
0.447 2471E-06 


0 . 5002033E-04 
0.1728000E-06 
0. 3252706E-06 
0.4777413E-06 


0.5082351 E-07 
0.2032941E-06 
0. 3557647E-06 
0. 508235  IE-06 


0.8131758E-07 
0.2337882E-06 
0.3862589E-06 
0. 5387294E-06 


0 . 1 1 18 l 1 7E-06 
0.2642823E-06 
0.4167530E-06 
0.5692230E-06 


CD 

<J1 


0.5997176E-06 
ROW  3 

0.0000000E  00 
0. 528 5647 E- 06 
0. 1 138446E-05 
0 • 1 748329E-05 
0.235821  IE-05 

ROW  4 

O.OOOOOOOE  00 
0. 1097787E-05 
0.2470023E-05 
0.3842259E-05 
0.5214493E-05 

ROW  5 

O.OOOOOOOE  00 
0 . 1 788988E-05 
0.4228518E-05 
0.6668048E-05 
0.9 107576E-05 

ROW  6 

O.OOOOOOOE  00 
0.5254U7E-04 
0.6352941E-05 
0.101 6470E-04 
0.1397646E-04 


0.63021 12E-06 


0 . 5082351E-07 
0.6505411E-06 
0.1260422E-05 
0. 1870305E-05 
0.2480186E-05 


0.8  131758E-07 
0. 1372235E-05 
0. 2 7 44470 E -05 
0.4116706E-05 
0. 5488939E-05 


0. 1 118117E-06 
0 . 2276894E-05 
0 . 47 16424E-05 
0.7155954E-05 
0.9595473E-05 


0. 1423058E-06 
0.3303528E-05 
0.71 15294E-05 
0.1092706E-04 
0. 1473880E-04 


0.6607058E-06 


0.5016263E-04 
0.7725175E-06 
0. 1382399E-05 
0. 1992281E-05 
0.2602163E-05 


0.28461 18E-06 
0.1646682E-05 
0.3018918E-05 
0. 4391 150E-05 
0 . 5763387E-05 


0.4065882E-06 
0. 2764800E-05 
0. 5204330E-05 
0.7643855E-05 
0. 1008339E-04 


0 . 5285647E-06 
0.4065882E-05 
0 . 7877644E-05 
0.1 1 68940E-04 
0. 15501 1 5E-04 


0.6911994E-06 


0.28461 18E-06 
0.8944938E-06 
0. 1504375E-05 
0.2114258E-05 
0.2724139E-05 


0 . 5054889E-04 
0.1921129E-05 
0 . 3293365E-05 
0.4665600E-05 
0.6037832E-05 


0.8233413E-06 
0. 3252706E-05 
0.5692236E-05 
0.8131765E-05 
0. 1057128E-04 


0. 1097787E-05 
0 . 4828235E-05 
0.86399 97 E-05 
0. 1245176E-04 
0. 1626350E-04 


0.72 16940 E-06 


0.4065882E-06 
0. 1016469E-05 
0.1626352E-05 
0. 2236233E-05 
0.2846 1 16E-05 


0.8233413E-06 
0.2195577E-05 
0.3567811E-05 
0.4940044 E-05 
0.6312281 E-05 


0.5130107  E-04 
0.3740612  E-05 
0.6180 142E-05 
0.8619662E-05 
0.1 105920E-04 


0.1788988E-05 
0. 5590588  E-05 
0.9402351 E-05 
0.  1321410E-04 
0. 1702584E-04 


¥ 


ROW  7 


O.OOOOOOOE  00 
0.3303528E-05 
0.8782305E-05 
0.1427124E-04 
0.1976016E-04 


0. 1728000E-06 
0.5439114E-04 
0 . 9880087E-05 
0.1536903E-04 
0.2085796E-04 


0.6505411E-06 
0.5488941E-05 
0. 1097788E-04 
0.1646680E-04 
0.2195574E-04 


0.  1372235E-05 
0. 6586730E-05 
0. 1207567E-04 
0.1756460E-04 
0. 2305352E-04 


0.2276894E-05 
0.7684514E-05 
0.1317344E-04 
0.1866238E-04 
0.2415 132E-04 


ROW  8 


O.OOOOOOOE  00 
0.4065882E-05 
0.1145562E-04 
0.1892669E-04 
0.2639773E-04 


0.2032941E-06 
0 .548894 1 E-05 
0.1294982E-04 
0.2042089E-04 
0. 2789193E-04 


0.7725175E-06 
0.5697299E-04 
0. 1444404E-04 
0.2191508E-04 
0. 29386 15E-04 


0. 16466 82 E-05 
0.8467 l 99C -05 
0. 1593825E-04 
0.2340930E-04 
0. 3088034E-04 


0.2764800E-05 
0.9961 411 E-05 
0. 1743247E-04 
0.2490351 E-04 
0.3237458E-04 


ROW  9 

O.OOOOOOOE  00 
0.4828235E-05 
0 . 1 43 1 190E-04 
0.2407002E-04 
0.33828 12E-04 


0. 2337882E-06 
0.6586730E-05 
0. 1626351E-04 
0.2602163E-04 
0. 3577973E-04 


0.8944938E-06 
0 • 8467 199E-05 
0. 1821514E-04 
0 . 2797325E-04 
0.3773137E-04 


0. 1921 1 29E-05 
0.6040865E-04 
0. 20 16 676 E-04 
0.2992488E-04 
0.  3 9 68  2 97  E-04 


0.3252706E-05 
0.1236027E-04 
0.2211839E-04 
0.3 187 648 E-04 
0 .4 163462 E-04 


ROW  10 

O.OOOOOOOE  00 
0. 5590588E-05 
0. 1729015E-04 
0.2964027E-04 
0.4199038E-04 

ROW  11 


0.2642823E-06 
0 .76845 14E-05 
0. 1976016E-04 
0.321 1029E-04 
0.4446038E-04 


0.1016469E-05 
0 .996141  IE-05 
0 . 2223020E-04 
0 • 3458030E-04 
0.4 6 93 04 2 E-04 


0. 2195577E-05 
0. 1236027E-04 
0.2470022E-04 
0. 3705033E-04 
0.494 0042 E-04 


0. 3740612E-05 
0.6482012E-04 
0.2717024E-04 
0.3952034E-04 
0. 5 187 046 E-04 


co 

-j 


O.OOOOOOOE  oc 
0.6352941E-05 
0.7032939E-04 
0.355764 6E-04 
0.5082351E-04 


0.2947765E-06 
0.8782305E-05 
0.2337880E-04 
0. 3862588E-04 
0. 5387289E-04 


0. 1138446E-05 
0.1 145562E-04 
0.2642823E-04 
0.4167526E-04 
0.5692233E-04 


0.2470023E-05 

0.1431190E-04 

0.2947764E-04 

0.4472470E-04 

0.5997172E-04 


0.4228518E-05 
0.1729015E-04 
0.3252705E-04 
0.477 7408 E-04 
0.63021 14E-04 


ROW  12 


O.OOOOOOOE  00 
0.7115294E-05 
0.2337880E-04 
0.418  1759E-04 
0.6026651E-04 

ROW  13 

O.OOOOOOOE  00 
0.7877644E-05 
0.2 64 28 2 3 E-04 
0.48302 68 E-04 
0.702  5844E-04 


0.3252706E-06 

0.9880087E-05 

0.7705843E-04 

0.4550737E-04 

0.6395628E-04 


0.3557647E-06 

0.1097788E-04 

0.3074821E-04 

0.5269384E-04 

0.7464956E-04 


0.1260422E-05 
0. 1294982E-04 
0.3074821E-04 
0.4919713E-04 
0.67 64609 E-04 


0. 1382399E-05 
0. 1444404E-04 
0.8512923E-04 
0. 570849 5E-04 
0.7 9 04074 E-04 


0.2744470E-05 
0. 1626351E-04 
0. 3443801 E-04 
0.5288694G-04 
0. 71 33585E-04 


0. 301891 8E-05 
0. 1821514E-04 
0. 3952038E-04 
0. 61 476 14 E-04 
0.8343186E-04 


0.4716424E-05 
0.1976016E-04 
0 .38 12779  E-04 
0.5657  670E-04 
0.7 502 566 E-04 


0. 5204330E-05 
0.2223020E-04 
0.4391 152E-04 
0.6586725E-04 
0. 8782303E-04 


ROW  14 

O.OOOOOOOE  00 
0.863 99 97E-05 
0.294 7764 E-04 
0. 5497073E-04 
0.8073825E-04 

ROW  15 

O.OOOOOOOE  00 
0.9402351E-05 


0.3862589E-06 
0. 1207567E-04 
0.3443801E-04 
0.60 12425E-04 
0.8589170E-04 


0.4167530E-06 
0. 1317344E-04 


0.1 504375E-05 
0. 1593825E-04 
0 . 3952038E-04 
0. 652777  IE-04 
0.9104525E-C4 


0. 1626352E-05 
0. 1743247E-04 


0. 3293365E-05 
0.2016676E-04 
0.9466373E-04 
0. 7043124E-04 
0.9619871E-04 


0. 356781 IE-05 
0.22 1 1839E-04 


0. 5692236E-05 
0.2470022E-04 
0.4981723E-04 
0.755847 1 E-04 
0.1013523E-03 


0.6180142E-05 
0. 2717024E-04 


oo 

oo 


0.3252705E-04 
0.617 6077E-04 
0. 9 16 4499 E- 04 

ROW  16 

O.OOOOOOOE  00 
0. 1 01 6470E-04 
0. 3557646E-04 
0. 1 1861 18E-03 
0. 1029176E-03 

ROW  17 

O.OOOOOOOE  00 
0.1092706E-04 
0.3862588E-04 
0.7547297E-04 
0.1144952E-03 

ROW  18 

O.OOOOOOOE  00 
0.1168940E-04 
0.4 167526E-04 
0.8233408E-04 
0.1 263166E-03 

ROW  19 

O.OOOOOOOE  00 
0.1 245176E-04 
0.4472470E-04 
0.8919531E-04 


0 . 38 12779E-04 
0.67737 62 E- 04 
0.9762179E-04 


0.4472 47  IE- 06 
0. 1427 124E-04 
0.4181759E-04 
0.7547297E-04 
0 . 1097788E-03 


0.4777413E-06 
0.1 536903E-04 
0.4550737E-04 
0.1 332693E-03 
0. 1223016E-03 


0 . 508235  IE-06 
0.1646680E-04 
0.4919713E-04 
0.9107571E-04 
0. 1351293E-03 


‘0.5387294E-06 
0. 1756460E-04 
0.5288694E-04 
0.9888226E-04 


0.4391152E-04 
0 . 7371441E-04 
0.1035987E-03 


0 . 1748329E-05 
0.1892669E-04 
0.4830268E-04 
0. 8233408E-04 
0 • 1 166400E-03 


0.  1870305E-05 
0.2042089E-04 
0.5269384E-04 
0.9107571E-04 
0. 1301082E-03 


0 . 199228  IE-05 
0.2191508E-04 
0 . 5708495E-04 
0 • 1498782E-03 
0 • 1439421E-03 


0.21 14258E-05 
0.2340930E-04 
0.6147614E-04 
0. 1086910E-03 


0. 4981723E-04 
0.7969131E-04 
0.1095755E-03 


0.3842259E-05 
0.2407002E-04 
0. 5497 073 E -04 
0.8919531E-04 
0. 1235010E-03 


0.4116706E-05 
0.2602163E-04 
0.6012425E-04 
0.9888226E-04 
0. 1379146E-03 


0.4391 150E-05 
0.2797325E-04 
0.6527771E-04 
0. 1086910E-03 
0. 1527549E-03 


0.4665600E-05 

0.2992488E-04 

0.7043124E-04 

0.1685610E-03 


0.1057839E-03 
0.8566810E-04 
0. 1 1 55524E-03 

4 


0.6668048  E-05 
0.2964027E-04 
0. 6176077E-04 
0. 9605 644 E -04 
0.1303622E-03 


0. 7 155 954 E-05 
0.321 1029E-04 
0.6773762  E-04 
0. 1066886E-03 
0.1457212E-03 


0 . 7643855  E-05 
0. 3458030E-04 
0 .737 1441 E-04 
0.117  5037  E-03 
0. 1615677E-03 


0.8131765E-05 
0. 3705033E-04 
0.7969131E-04 
0. 1284410E-03 


0.1383211E-03 
ROW  20 

O.OOOOOOOE  00 
0.132 1410E-04 
0.4777408E-04 
0.9605644E-04 
0. 1 504476E-03 

ROW  21 

O.OOOOOOOE  00 
0.1397646E-04 
0.5082351E-04 
0 . 1029 176E-03 
0.2 126351E-03 

ROW  22 

O.OOOOOOOE  00 
0.147  3880E-04 
0.5387289E-04 
0.1097788E-03 
0. 1 748324E-03 

ROW  23 

O.OOOOOOOE  00 
0.15501 1 5E-04 
0.5692233E-04 
0.1 166400E-03 
0.1 870302E-03 


0.1482012E-03 


0.5692230E-06 
0.1866238E-04 
0.5657 670E-04 
0.1066886E-03 
0.1614558E-03 


0. 5997176E-06 
0.1976016E-04 
0.6026651E-04 
0 . 1 144952E-03 
0. 1748324E-03 


0 . 6302 112E-06 
0. 2085796E-04 
0.6395628E-04 
0.1223016E-03 
0. 2382702E-03 


0.6607058E-06 
0.2195574E-04 
0.6764609E-04 
0.1301082E-03 
0.2017 181E-03 


0. 1580813E-03 


0.2236233E-05 
0.2490351E-04 
0.6586725E-04 
0. 1175037E-03 
0 . 1 724642E-03 


0.23582 1 IE-05 
0.2639773E-04 
0.7025844E-04 
0. 1263166E-03 
0. 1870302E-03 


0. 2480186E-05 
0.2789193E-04 
0.7464956E-04 
0. 1351293E-03 
0 . 20 17 18  IE-03 


0.2602163E-05 
0 . 2938615E-04 
0.7904074E-04 
0. 1439421E-03 
0. 266467  2E-03 


0. 1679613E-03 


0.4940044E-05 
0. 3187648E-04 
0.7558471E-04 
0.1284410E-03 
0.1834724E-03 


0. 52 14493 E -05 
0.3382812E-04 
0.8073825E-04 
0. 1383211E-03 
0.1992277E-03 


0.5488939E-05 
0.3577973E-04 
0.8589170E-04 
0. 1482012E-03 
0.2151659E-03 


0. 5763387E-05 
0. 3773 137E-04 
0.9104525E-04 
0.1580813E-03 
0. 2312262E-03 


0. 1778415E-03 


0.8619662E-05 
0.3952034E-04 
0.8566810E-04 
0. 1894390E-03 
0.1944808E-03 


0.9107  576E-05 
0 .4199038  E-04 
0.9164499E-04 
0. 1504476E-03 
0.2114255E-03 


0.9595473E-05 
0.4446038  E-04 
0.9762179E-04 
0. 1614558E-03 
0.2286140E-03 


0. 1008339E-04 
0.4693042E-04 
0 . 1035987  E-03 
0. 1724642E-03 
0. 2459851 E-03 


ID 
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ROW  24 


O.OOOOOOOE  00 
0. 1626350E-04 
0.5997172E-04 
0.12350 10E-03 
0.1992277E-03 

ROW  25 

O.OOOOOOOE  00 
0.1702584E-04 
0.63021 14E-04 
0. 1 303622E-03 
0.2114255E-03 


0.6911994E-06 
0.2305352E-04 
0. 7 133585E-04 
0.1 379 146E-03 
0.2  151659E-03 


0. 72 16940E-06 
0.2415132E-04 
0. 75025 66E- 04 
0.1457212E-03 
0.2286140E-03 


0. 2724139E-05 
0.3088034E-04 
0.8343186E-04 
0.1 527549E-03 
0.2312262E-03 


0.2846 1 16E-05 
0 • 3237458E-04 
0. 8782303E-04 
0. 1615677E-03 
0.2459851E-03 


0. 6037 832 E -05 
0.39 68297 E-04 
0.9619871E-04 
0. 1679613E-03 
0. 2973471 E-03 


0. 6312281E-05 
0 . 41 634 62 E-04 
0. 1013523E-03 
0.1778415E-03 
0.2634786E-03 


0.1057128E-04 
0.4940042 E-04 
0. 1095755E-03 
0.1834724E-03 
0.2634786E-03 


0. 1 105920E-04 
0.5187046E-04 
0.1 155524E-03 
0.1944808E-03 
0.3310333E-03 


POINTS  IN  CONTACT  OBTAINEO  FROM  QUADRATIC  PROGRAMMING. 
12345678  9101112 
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